Modelos de regresión para series de tiempo

Prerrequisitos

Se requieren cargar las siguientes paqueterías para seguir los ejemplos.

library(easypackages)
libraries("tidyverse","fpp3","plotly", "patchwork")

Introducción

Los modelos que veremos aquí tienen como idea principal encontrar relaciones lineales entre la serie que queremos pronosticar, \(y\), con una o más series distintas, x. En otras palabras, pronosticaremos los valores futuros de una serie, a partir de los cambios en otra serie que la afecte.

Es muy común querer predecir de esta forma. Por ejemplo, una tienda de helados podría encontrar una relación entre sus ventas (\(y\)) y la temperatura (\(x_1\)). O las ventas de Nike, a partir de cuánto gastan en publicidad y mercadotecnia.

En la literatura podemos encontrar muchos nombres para las variables \(y\) ^ \(x\). P. ej.

\(y\) (var. de pronóstico)	\(x\) (vars. predictoras)
Var. dependiente	Vars. independientes
Explicada	Explicativas
Regresada	Regresoras
Respuesta	Estímulo
Resultado	Covariante
Controlada	De control

El modelo lineal

El caso más sencillo sería un modelo de regresión lineal simple, de la forma:

\[ y_t = \beta_0 + \beta_1 x_t + \varepsilon_t \]

donde

• \(\beta_0\) es conocido como el intercepto y representa el valor predicho cuando \(x = 0\).

• \(\beta_1\) es la pendiente de la recta. Nos indica el cambio promedio en \(y\), ante un cambio en una unidad de \(x\).

• El término de error, \(\varepsilon_t\) se asume aleatorio y decimos que captura los cambios debido a todas las otras variables que pudieran llegar a afectar a \(y_t\), que no están explícitamente especificadas en el modelo.

La recta resultante está dada entonces por \(\beta_0 + \beta_1 x_t\), y la diferencia que existe en los puntos reales y ésta es \(\varepsilon_t\).

(Hyndman, 2019)

Ejemplo: gasto de consumo en EEUU

Como primer ejemplo, veamos las tasas de crecimiento del gasto de consumo, \(y\), y su relación con el ingreso personal disponible, \(x\).

La gráfica de tiempo de ambas series:

us_change %>%
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Consumption, colour = "Consumo")) +
  geom_line(aes(y = Income, colour = "Ingreso")) +
  ylab("cambio %") + xlab("Año") +
  guides(colour=guide_legend(title="Series")) + 
  theme(legend.position = "top")

Un diagrama de dispersión entre ambas series, para ver una posible correlación.

us_change %>%
  ggplot(aes(x=Income, y=Consumption)) +
    ylab("Consumo (cambio % trimestral)") +
    xlab("Ingreso (cambio % trimestral)") +
    geom_point() +
    geom_smooth(method="lm", se=FALSE)

La línea azul en el gráfico sigue la ecuación que describe la regresión lineal que tiene por variable dependiente al consumo y como independiente al ingreso.

¿Qué es una regresión lineal?

El término de regresión fue acuñado por primera vez por Francis Galton en 1886. Él estaba estudiando la relación que existe entre la estatura de los hijos y la estatura de los padres.

Lo que encontró fue lo siguiente, en resumen:

• Los padres más altos, tendían a tener hijos más altos, mientras que los padres bajos tendían a tener hijos bajos.

• En promedio, los hijos de padres altos no logran ser más altos que ellos. Similarmente, los hijos de padres bajos, en promedio son más altos que sus papás.

• Así, Galton decía que había una tendencia a regresar a la estatura promedio.

De no cumplirse la regresión de Galton, sería común tener gente de la estatura de un Hobbit, y también de la estatura de un gigante.

Entonces, el análisis de regresión en tiempos modernos trata sobre la relación de la dependencia entre una variable \(y\), respecto de una o más variables exógenas (regresoras \(x\)) para predecir el valor promedio de la variable dependiente.

Regresión y causalidad

“Una relación estadística, por más fuerte y sugerente que sea, nunca podrá establecer una conexión causal: nuestras ideas de causalidad deben provenir de estadísticas externas y, en último término, de una u otra teoría” (Kendall & Stuart, 1961)

Regresión \(\neq\) Causalidad

Correlación \(\neq\) Causalidad

• Una relación estadística por sí misma no puede implicar causalidad.
  • Se debe acudir a consideraciones a priori o teóricas.

• La causalidad puede determinarse también por sentido común.

¿Qué significa que un modelo sea lineal?

Se puede hablar de linealidad en dos sentidos:

1. Linealidad en las variables; \(x\).

• La esperanza condicional de \(y\) es una función lineal de \(x_i\).

2. Linealidad en los parámetros; \(\beta\).

• La esperanza condicional de \(y\) es lineal en los parámetros \(\beta_i\).

Para un modelo de regresión lineal, solo nos interesa que sea lineal en los parámetros.

Todas estas funciones son lineales en los parámetros y pueden ser estimadas mediante un modelo de regresión lineal

Así, un modelo de regresión lineal puede generar una recta, o una variedad de curvas, dependiendo la forma funcional que se elija.

Mínimos cuadrados ordinarios (MCO-OLS)

\[ H_0: \beta_0 = 0 \\ H_1: \beta_0 \neq 0 \] \[ H_0: \beta_1 = 0 \\ H_1: \beta_1 \neq 0 \]

\[ \hat{y}_t = \hat{\beta}_0 + \hat{\beta}_1x_t + \hat{\varepsilon}_t \\ \hat{y}_t = 0.54454 + 0.27183 x_t + \hat{\varepsilon}_t \]

\[ y_{consumo} = \beta_0 + \beta_1 x_{income} + \beta_2 x_{production} + \beta_3 x_{savings} + \beta_4 x_{unemployment} \]

Ejemplos

US % change

us_change

us_change %>% 
  as_tibble() %>% 
  select(-Quarter) %>% 
  GGally::ggpairs()

Una correlación, por más fuerte que sea entre dos variables, no puede implicar por sí misma causalidad.

us_change %>% 
  pivot_longer(cols = -Quarter) %>% 
  ggplot(aes(x = Quarter, y = value, color = name)) +
  geom_line() +
  facet_wrap(~ name, scales = "free_y") +
  theme(legend.position = "none")

Graficando nuestra variable de pronóstico (Consumo) vs. cada una de las variables predictoras:

us_change %>% 
  pivot_longer(cols = -c(Quarter, Consumption)) %>% 
  ggplot(aes(x = Quarter, y = value, color = name)) +
  geom_line() +
  geom_line(aes(y = Consumption), color = "black") +
  facet_wrap(~ name, scales = "free_y") +
  theme(legend.position = "none")

Regresión lineal simple

Realizamos un primer modelo, donde utilizaremos de variable predictora al ingreso disponible, para pronosticar el consumo.

\[ y_{t,Consumo} = \beta_0 + \beta_{1}x_{t,Ingreso} + \varepsilon_t \]

fit1 <- us_change %>% 
  model(reg_lin_simple = TSLM(Consumption ~ Income)
        )
fit1 %>%  report()

Series: Consumption 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-2.58236 -0.27777  0.01862  0.32330  1.42229 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.54454    0.05403  10.079  < 2e-16 ***
Income       0.27183    0.04673   5.817  2.4e-08 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5905 on 196 degrees of freedom
Multiple R-squared: 0.1472, Adjusted R-squared: 0.1429
F-statistic: 33.84 on 1 and 196 DF, p-value: 2.4022e-08

Prueba de significancia individual para cada uno de los parámetros (\(\beta_i\))

\[ H_0: \beta_i = 0 \]

Prueba de significancia conjunta (para revisar si nuestro modelo sirve o no)

\[ H_0: \beta_1 = \beta_2 = \beta_3 = \ldots = 0 \]

augment(fit1) %>% 
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Consumption, color = "Datos")) +
  geom_line(aes(y = .fitted, color = "Fitted"))+
  xlab("Año") + ylab(NULL) +
  ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
  guides(color = guide_legend(title = NULL))

El modelo no parece capturar adecuadamente la variación de los datos reales.

augment(fit1) %>% 
  ggplot(aes(x = Consumption, y = .fitted)) +
  geom_point() +
  ylab("Fitted (valores ajustados)") +
  xlab("Datos (reales históricos)") +
  ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
  geom_abline(intercept = 0, slope = 1)

fit1 %>% 
  gg_tsresiduals()

augment(fit1) %>% 
  features(.resid, ljung_box, lag= 10, dof = 2)

Es evidente que este modelo se puede mejorar. Probemos incluyendo las otras predictoras.

Regresión lineal múltiple

Podríamos proponer ahora un modelo en donde el consumo sea una función del ingreso, la producción, los ahorros y el desempleo:

\[ y_{t,Consumo} = \beta_0 + \beta_{1}x_{t,Ingreso} + \beta_2x_{t,Producción} + \beta_3x_{t,Ahorros} + \beta_4x_{t,Desempleo} + \varepsilon_t \]

fit2 <- us_change %>% 
  model(
    reg_lin_multiple = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
  )
report(fit2)

Series: Consumption 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-0.90555 -0.15821 -0.03608  0.13618  1.15471 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.253105   0.034470   7.343 5.71e-12 ***
Income        0.740583   0.040115  18.461  < 2e-16 ***
Production    0.047173   0.023142   2.038   0.0429 *  
Savings      -0.052890   0.002924 -18.088  < 2e-16 ***
Unemployment -0.174685   0.095511  -1.829   0.0689 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3102 on 193 degrees of freedom
Multiple R-squared: 0.7683, Adjusted R-squared: 0.7635
F-statistic:   160 on 4 and 193 DF, p-value: < 2.22e-16

augment(fit2) %>% 
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Consumption, color = "Datos")) +
  geom_line(aes(y = .fitted, color = "Fitted"))+
  xlab("Año") + ylab(NULL) +
  ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
  guides(color = guide_legend(title = NULL))

Este modelo parece capturar más variación de los datos históricos.

augment(fit2) %>% 
  ggplot(aes(x = Consumption, y = .fitted)) +
  geom_point() +
  ylab("Fitted (valores ajustados)") +
  xlab("Datos (reales históricos)") +
  ggtitle("Cambios porcentuales en el gasto de Consumo en EEUU") +
  geom_abline(intercept = 0, slope = 1)

fit2 %>% 
  gg_tsresiduals()

augment(fit2) %>% 
  features(.resid, ljung_box, lag= 10, dof = 2)

df <- left_join(us_change, residuals(fit2), by = "Quarter")
df %>% 
  select(-c(Consumption, .model)) %>% 
  pivot_longer(cols = c(Income:Unemployment)) %>% 
  ggplot(aes( x = value, y = .resid, color = name)) + 
  geom_point() + ylab("Residuales") + xlab("Predictoras") +
  facet_wrap(~ name, scales = "free_x") +
  theme(legend.position = "none")

augment(fit2) %>% 
  ggplot(aes(x = .fitted, y = .resid)) +
  geom_point() +
  labs(x = "Ajustados", y = "Residuales")

glance(fit2) %>% 
  select(adj_r_squared, AIC, AICc, BIC)

Selección de predictoras

1. Escoger subconjuntos de predictoras y probarlo.

fit3 <- us_change %>% 
  model(r1 = TSLM(Consumption ~ Income),
        r2 = TSLM(Consumption ~ Income + Production),
        r3 = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
        r4 = TSLM(Consumption ~ Income + Production + Savings),
        r5 = TSLM(Consumption ~ Income + Savings + Unemployment),
        r6 = TSLM(Consumption ~ Income + Production + Unemployment),
        r7 = TSLM(Consumption ~ Income + Savings)
        )
fit3 %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

fit3 %>% 
  select(r3) %>% 
  report()

Series: Consumption 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-0.90555 -0.15821 -0.03608  0.13618  1.15471 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.253105   0.034470   7.343 5.71e-12 ***
Income        0.740583   0.040115  18.461  < 2e-16 ***
Production    0.047173   0.023142   2.038   0.0429 *  
Savings      -0.052890   0.002924 -18.088  < 2e-16 ***
Unemployment -0.174685   0.095511  -1.829   0.0689 .  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3102 on 193 degrees of freedom
Multiple R-squared: 0.7683, Adjusted R-squared: 0.7635
F-statistic:   160 on 4 and 193 DF, p-value: < 2.22e-16

2. Backwards stepwise regression:

• Empezamos con un modelo que contenga todas las predictoras.
• Quitamos una a la vez.
• Mantenemos el modelo si mejora la medida de desempeño predictivo (\(\bar{R}^2\), AICc,…).
• Seguirlo haciendo hasta no encontrar mejoras adicionales.

Selección del modelo a través de backwards stepwise regression con base en el \(AICc\) sugiere conservar el modelo que incluye todas las predictoras.

us_change %>% 
  model(i = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
        ii = TSLM(Consumption ~ Income + Production + Savings),
        iii = TSLM(Consumption ~ Income + Production + Unemployment),
        iv = TSLM(Consumption ~ Income + Savings + Unemployment),
        v = TSLM(Consumption ~ Production + Savings + Unemployment)
        ) %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

Con base en \(BIC\) parece que el mejor modelo es el que incluye todas las predictoras menos el desempleo.

us_change %>% 
  model(i = TSLM(Consumption ~ Income + Production + Savings + Unemployment),
        ii = TSLM(Consumption ~ Income + Production + Savings),
        iii = TSLM(Consumption ~ Income + Production),
        iv = TSLM(Consumption ~ Income + Savings),
        v = TSLM(Consumption ~ Production + Savings),
        vi = TSLM(Consumption ~ Income + Savings + Unemployment),
        vii = TSLM(Consumption ~ Production + Savings + Unemployment),
        viii = TSLM(Consumption ~ Income + Production + Unemployment)
        ) %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

3. Forwards stepwise regression:

• Comenzar con un modelo que solo incluya al intercepto.
• Se van agregando las predictoras una a la vez.
• La predictora que mejore más al modelo se mantiene.
• Se itera hasta no tener mejoría adicional.

us_change %>% 
  model(i = TSLM(Consumption ~ Income),
        ii = TSLM(Consumption ~ Production),
        iii = TSLM(Consumption ~ Savings),
        iv = TSLM(Consumption ~ Unemployment)
        ) %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

El mejor modelo hasta ahora es el que toma de predictora a la producción.

us_change %>% 
  model(
    o   = TSLM(Consumption ~ Production),
    i   = TSLM(Consumption ~ Production + Income),
    ii  = TSLM(Consumption ~ Production + Savings),
    iii = TSLM(Consumption ~ Production + Unemployment)
  ) %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

El mejor hasta ahora es el modelo con la producción y el ingreso.

us_change %>% 
  model(
    i   = TSLM(Consumption ~ Production + Income),
    ii  = TSLM(Consumption ~ Production + Income + Savings),
    iii = TSLM(Consumption ~ Production + Income + Unemployment)
  ) %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

El modelo mejora agregando los ahorros también.

us_change %>% 
  model(
    i   = TSLM(Consumption ~ Production + Income + Savings),
    ii  = TSLM(Consumption ~ Production + Income + Savings + Unemployment)
  ) %>% 
  glance() %>% 
  select(.model, adj_r_squared, AIC, AICc, BIC)

Utilizando la adj_r_squared o el AICc, concluimos lo mismo que antes: el modelo que incluye a todas las variables es el mejor. ### {-}

Pronóstico

Pronósticos ex-ante

En estos pronósticos solo se utiliza información disponible hasta el último dato del histórico. A estos pronósticos se les considera como pronósticos reales. Aquí las predictoras se deben pronosticar antes de poder producir el pronóstico de la variable de interés.

Pronósticos ex-post

Con estos pronósticos se utiliza información real disponible de las predictoras. Estos pronósticos ya no son reales (en el sentido estricto). La variable a pronosticar (\(y\)) sigue siendo desconocida.

Pronósticos basados en escenarios

fit_escenarios <- us_change %>% 
  model(lineal = TSLM(Consumption ~ Income + Savings + Unemployment))

# Necesitamos agregar nuevos datos de las predictoras
escenarios <- scenarios(
  optimista = new_data(us_change, 4) %>% 
    mutate(Income = c(0,0.2,0,1.2), Savings = c(0,-0.1,0.5,-1), Unemployment = -0.1),
  pesimista = new_data(us_change, 4) %>% 
    mutate(Income = -1, Savings = -0.5, Unemployment = 0.1),
  names_to = "Escenario"
)

fc_escenarios <- fit_escenarios %>% 
  forecast(new_data = escenarios)

us_change %>% 
  autoplot(Consumption) +
  autolayer(fc_escenarios)

Vamos a generar un pronóstico ex-ante del consumo. Para esto, necesitamos primero producir pronósticos de las predictoras:

# mod_predictoras <- function(predictora, horizonte = 4) {
#   us_change %>% 
#     model(predictora = ARIMA(as.formula(predictora)) %>% 
#     forecast(h = horizonte)
# }
# 
# mod_predictoras(predictora = Income)

ingreso <-  us_change %>% 
  model(ETS = ETS(Income),
        ARIMA = ARIMA(Income)
        ) %>% 
  forecast(h = 4) 

ingreso %>% 
  autoplot(us_change, level = NULL)

produccion <- us_change %>% 
  model(
    ETS = ETS(Production),
    ARIMA = ARIMA(Production)
  ) %>% 
  forecast(h = 4)
produccion %>% 
  autoplot(us_change, level = NULL)

ahorro <- us_change %>% 
  model(
    ETS = ETS(Savings),
    ARIMA = ARIMA(Savings)
  ) %>% 
  forecast(h = 4)
ahorro %>% 
  autoplot(us_change, level = NULL)

desempleo <- us_change %>% 
  model(
    ETS = ETS(Unemployment),
    ARIMA = ARIMA(Unemployment)
  ) %>% 
  forecast(h = 4)
desempleo %>% 
  autoplot(us_change, level = NULL)

Teniendo ya los pronósticos de cada predictora, podemos proceder a generar el pronóstico del consumo.

fit <- us_change %>% 
  model(
    `Regresión lineal múltiple` = TSLM(Consumption ~ Income + Production + Savings + Unemployment)
  )

datos_futuros <- new_data(us_change,4) %>% 
  mutate(Income = ingreso %>% filter(.model == "ARIMA") %>% pull(.mean), 
         Savings = ahorro %>% filter(.model == "ARIMA") %>% pull(.mean), 
         Unemployment = desempleo %>% filter(.model == "ARIMA") %>% pull(.mean),
         Production = produccion %>% filter(.model == "ARIMA") %>% pull(.mean))

datos_futuros

fc <- forecast(fit, datos_futuros)

fc %>% 
  autoplot(us_change)

fc %>% 
  autoplot(us_change %>% filter_index("2016 Q1" ~ .))

Inclusión de predictoras útiles: Producción de cerveza

recent_production <- aus_production %>% 
  filter(year(Quarter) >= 1992)

recent_production

recent_production %>% 
  autoplot(Beer) +
  labs(x = "Año", y = "Megalitros", 
       title = "Producción de cerveza trimestral en Australia")

Existen varias predictoras que pueden ser útiles en el análisis de regresión.

1. Predictora de tendencia trend()

\[ y_t = \beta_0 + \beta_1t + \varepsilon_t \]

2. Variables dummy estacionales

Las variables dummy (también llamadas dicotómicas, binarias, …) toman solo dos valores: 1 o 0 (verdadero/falso).

Para agregar variables dummy estacionales, basta con escribir season().

recent_production %>% 
  select(Quarter,Beer) %>% 
  mutate(tendencia = seq_along(recent_production$Quarter),
         q2 = if_else(quarter(Quarter)==2,1,0),
         q3 = if_else(quarter(Quarter)==3,1,0),
         q4 = if_else(quarter(Quarter)==4,1,0)
         ) %>% 
  model(TSLM(Beer ~ tendencia + q2 + q3 + q4)) %>% 
  report()

Series: Beer 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-42.9029  -7.5995  -0.4594   7.9908  21.7895 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 441.80044    3.73353 118.333  < 2e-16 ***
tendencia    -0.34027    0.06657  -5.111 2.73e-06 ***
q2          -34.65973    3.96832  -8.734 9.10e-13 ***
q3          -17.82164    4.02249  -4.430 3.45e-05 ***
q4           72.79641    4.02305  18.095  < 2e-16 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.23 on 69 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9199
F-statistic: 210.7 on 4 and 69 DF, p-value: < 2.22e-16

fit_beer <- recent_production %>% 
  model(TSLM(Beer ~ trend() + season()))

report(fit_beer)

Series: Beer 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-42.9029  -7.5995  -0.4594   7.9908  21.7895 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)   441.80044    3.73353 118.333  < 2e-16
trend()        -0.34027    0.06657  -5.111 2.73e-06
season()year2 -34.65973    3.96832  -8.734 9.10e-13
season()year3 -17.82164    4.02249  -4.430 3.45e-05
season()year4  72.79641    4.02305  18.095  < 2e-16
                 
(Intercept)   ***
trend()       ***
season()year2 ***
season()year3 ***
season()year4 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.23 on 69 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9199
F-statistic: 210.7 on 4 and 69 DF, p-value: < 2.22e-16

p <- augment(fit_beer) %>% 
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Beer, color = "Datos")) +
  geom_line(aes(y = .fitted, color = "Fitted")) +
  labs(x = "Año", y = "Megalitros")

ggplotly(p)

Al ver la gráfica, vemos que existe un trimestre muy por debajo del resto. Podemos agregar una variable de intervención para ese periodo.

Spike variables

Capturan el efecto de un solo periodo.

cerveza <- recent_production %>% 
  select(Quarter, Beer) %>% 
  mutate(
    q2_94 = if_else(Quarter == yearquarter("1994 Q2"),1,0),
    q4_04 = if_else(Quarter == yearquarter("2004 Q4"),1,0)
  )
cerveza

Ajustamos un modelo, corrigiendo por ese periodo outlier (1994 Q2).

fit_beer <- cerveza %>% 
  model(TSLM(Beer ~ trend() + season() + q2_94 + q4_04))

report(fit_beer)

Series: Beer 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-27.0379  -5.8811  -0.6895   7.7209  21.7895 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)
(Intercept)   441.84123    3.32189 133.009  < 2e-16
trend()        -0.34137    0.05975  -5.713 2.78e-07
season()year2 -33.39369    3.55788  -9.386 7.83e-14
season()year3 -17.82164    3.55460  -5.014 4.16e-06
season()year4  75.32030    3.60790  20.876  < 2e-16
q2_94         -24.03384   11.24265  -2.138 0.036190
q4_04         -45.41027   11.15547  -4.071 0.000126
                 
(Intercept)   ***
trend()       ***
season()year2 ***
season()year3 ***
season()year4 ***
q2_94         *  
q4_04         ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 10.81 on 67 degrees of freedom
Multiple R-squared: 0.9426, Adjusted R-squared: 0.9375
F-statistic: 183.4 on 6 and 67 DF, p-value: < 2.22e-16

p <- augment(fit_beer) %>% 
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Beer, color = "Datos")) +
  geom_line(aes(y = .fitted, color = "Fitted")) +
  labs(x = "Año", y = "Megalitros")

ggplotly(p)

Cambios de nivel

Capturan el efecto a partir de cierto periodo.

cerveza <- cerveza %>% 
  mutate(d2000 = if_else(year(Quarter)>=2000,1,0))
cerveza

fit_beer <- cerveza %>% 
  model(TSLM(Beer ~ trend() + season() + d2000))

report(fit_beer)

Series: Beer 
Model: TSLM 

Residuals:
     Min       1Q   Median       3Q      Max 
-43.0235  -7.5945  -0.4539   7.7754  21.4829 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   441.9154     3.8644 114.354  < 2e-16 ***
trend()        -0.3548     0.1308  -2.712  0.00846 ** 
season()year2 -34.6452     3.9985  -8.665 1.36e-12 ***
season()year3 -17.8046     4.0536  -4.392 4.02e-05 ***
season()year4  72.8279     4.0594  17.941  < 2e-16 ***
d2000           0.7282     5.6402   0.129  0.89766    
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 12.32 on 68 degrees of freedom
Multiple R-squared: 0.9243, Adjusted R-squared: 0.9188
F-statistic: 166.1 on 5 and 68 DF, p-value: < 2.22e-16

p <- augment(fit_beer) %>% 
  ggplot(aes(x = Quarter)) +
  geom_line(aes(y = Beer, color = "Datos")) +
  geom_line(aes(y = .fitted, color = "Fitted")) +
  labs(x = "Año", y = "Megalitros")

ggplotly(p)

cerveza %>% 
  gg_tsdisplay(Beer %>% difference(4), plot_type = "partial")

tibble(mujer = c(1,0,1, 1),
       hombre = c(0,1,0, 0),
       nombre = c("Andrea","Juan","Sofía","Fer"))

Cuando se crean varaibles dummy, la cantidad de dummies a generar es una menos que el total de categorías. Si son dos categorías, solo necesitamos una dummy. Si son tres, necesitamos 2, etc.

Regresiones no lineales

Modelos exponenciales

• Modelos log-log

\[ \log{y_t} = \beta_0 + \beta_1 \log{x_t} \]

La interpretación de los coeficientes (\(\beta_s\)) es como elasticidades (cambios porcentuales).

• Modelos lin-log

\[ y_t = \beta_0 + \beta_1 \log{x_t} \]

• Modelos log-lin

\[ \log{y_t} = \beta_0 + \beta_1 x_t \]

Modelos de regresión lineal por partes (piecewise)

Aquí la regresión se calcula por partes en el tiempo, para capturar distintas tendencias.

boston_men <- boston_marathon %>% 
  filter(Event == "Men's open division") %>% 
  mutate(Minutes = as.numeric(Time)/60)

p <- boston_men %>% 
  autoplot(Minutes) + 
  ggtitle("Tiempos ganadores del maratón de Boston, categoría abierta de hombres")

ggplotly(p)

fit_boston <- boston_men %>% 
  model(
    lineal = TSLM(Minutes ~ trend()),
    )

fc_boston <- fit_boston %>% forecast(h = 10)

boston_men %>% 
  autoplot(Minutes) +
  geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
  autolayer(fc_boston, alpha = 0.5, level = 95) +
  ggtitle("Maratón de Boston, cat. abierta de hombres")

Vamos a modelar los tiempos con una regresión por partes.

fit_boston <- boston_men %>% 
  model(
    lineal = TSLM(Minutes ~ trend()),
    exponencial = TSLM(log(Minutes) ~ trend()),
    `Reg. por partes` = TSLM(Minutes ~ trend(knots = c(1940,1980)))
  )

fc_boston <- fit_boston %>% forecast(h = 10)

boston_men %>% 
  autoplot(Minutes) +
  geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
  autolayer(fc_boston, alpha = 0.5, level = 95) +
  ggtitle("Maratón de Boston, cat. abierta de hombres")

Intentamos ahora poniendo los cambios en la tendencia en otros periodos:

fit_boston <- boston_men %>% 
  model(
    lineal = TSLM(Minutes ~ trend()),
    exponencial = TSLM(log(Minutes) ~ trend()),
    `Reg. por partes` = TSLM(Minutes ~ trend(knots = c(1940,1980))),
    `Reg. por partes2` = TSLM(Minutes ~ trend(knots = c(1975))),
    `Reg. por partes3` = TSLM(Minutes ~ trend(knots = c(1915, 1950, 1988))),
    `Reg. por partes4` = TSLM(Minutes ~ trend(knots = c(1915, 1927, 1950, 1985)))
  )

fc_boston <- fit_boston %>% forecast(h = 10)

boston_men %>% 
  autoplot(Minutes) +
  geom_line(aes(y = .fitted, color = .model), data = fitted(fit_boston)) +
  autolayer(fc_boston, alpha = 0.5, level = 95) +
  ggtitle("Maratón de Boston, cat. abierta de hombres")

accuracy(fit_boston) %>% 
  arrange(MAPE)

Probamos transformando la serie con Box-Cox:

lambda <- boston_men %>% 
  features(Minutes, guerrero) %>% 
  pull(lambda_guerrero)

p1 <- boston_men %>% 
  autoplot(Minutes)

p2 <- boston_men %>% 
  autoplot(box_cox(Minutes, lambda = lambda))

p3 <- boston_men %>% 
  autoplot(log(Minutes))

p1/p2/p3